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Abstract 

It is demonstrated that the well-regularized hypergeometric functions can be evaluated 
directly and numerically. The package NumExp is presented for expanding hypergeomet- 
ric functions and/or other transcendental functions in a small regularization parameter. 
Hypergeometric function is expressed as a Laurent series in the regularization parameter 
and the coefficients are evaluated numerically by using multi-precision finite difference 
method. This elaborate expansion method works for a wide variety of hypergeometric 
functions, which are needed in the context of dimensional regularization for loop integrals. 
The divergent and finite parts can be extracted from the final result easily and simulta- 
neously. In addition, there is almost no restriction on the parameters of hypergeometric 
functions. 
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Program summary 

Program title: NumExp 
Catalogue identifier: 
Program summary URL: 
Program obtainable from: 
Licensing provisions: none 

No. of lines in distributed program, including test data, etc.: 219 
No. of bytes in distributed program, including test data, etc.: 9 071 
Distribution format: tar . gz 

Programming language: Mathcmatica and/or Python 

Computer: Any computer where Mathematica or Python is running. 

Operating system: Linux, Windows 

External routines: mpmath library (for Python) 

Classification: 4.4, 5, 11.1 

Nature of problem: Expansion of hypergeometric functions and/or other transcendental 
functions in a small parameter e. These expansions are needed in the context of 
dimensional regularization for loop integrals. 

Solution method: Hypergeometric function is expressed as a Laurent series in the regu- 
larization parameter e, where the coefficients are evaluated numerically by multi- 
precision finite difference method. 

Running time: Generally it is below a few seconds, depending on the complexity of the 
problem. 
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1. Introduction 



Hypcrgeometric functions and their extensions are used frequently in the calculation 
of Feynman integrals in quantum field theory [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. The hy- 
pergeometric representation makes the calculation of Feynman integrals less tricky and 
systematic. In the past decades, several methods, such as the Mellin-Barnes techniques 
[1, 11, 12], the negative dimensional integration method [13, 14, 5] or its optimized 
version, the method of brackets [15, 16, 17, 18], the differential equations techniques 
[19, 8], the DRA method [20, 21] and so on, have been used to obtain the hypergeometric 
representation of the Feynman integrals. The hypergeometric representation keeps the 
expression of the result in a compact form and spurious singularities cancel each other 
out automatically. Besides, the hypergeometric functions have good algebraic and ana- 
lytic properties, which make it easy to study the physical problems in different kinematic 
regions. 

In the context of dimensional regularization [22], the space-time dimension D = 
4 — 2e appears in the parameters of hypergeometric functions, where the parameter e 
regulates infrared and/or ultraviolet divergences. Formally, the solution in the form of 
hypergeometric functions can be expressed as a Laurent series in e, but in practice this 
e-expansion of hypergeometric functions is still not a trivial task. 

Recently, a number of elegant algorithms [23, 24] and packages like nestedsums 
[25, 26], Xsummer [27], HypExp [28, 29] as well as HYPERDIRE [30, 31] are developed 
to perform analytic e-expansion of hypergeometric or transcendental functions. These 
implementations, however, are restricted to the expansion of some special classes of hy- 
pergeometric functions about integer and/or half-integer parameters. It is increasingly 
obvious that more general and convenient algorithms are needed to meet the demand of 
practical calculations. 

In this work, a numerical algorithm is developed for the e-expansion of a wide variety 
of hypergeometric functions. Specifically, the generalized hypergeometric functions p F q , 
the Appell hypergeometric functions and the Horn-type hypergeometric functions of two 
variables, can all be expanded in e by this method. In principle hypergeometric func- 
tions are treated as mathematical objects and there is no restriction on the form of their 
parameters. Analytic continuation is performed automatically, which significantly sim- 
plifies the practical calculation. In addition, the regularized hypergeometric function can 
be evaluated numerically and lengthy expressions of the analytic expansion are avoided. 
Furthermore, no knowledge of harmonic or multiple polylogarithms [32, 33, 34, 35] or 
other newly defined special functions is needed. As a result, the well regularized hyper- 
geometric functions can be treated as common functions in the practical calculations. 

The numerical algorithm is based on the ansatz that the regularized hypergeometric 
function can be expressed as a Laurent series in the regularization parameter, where the 
coefficients are evaluated numerically by multi-precision finite difference method [36]. 
Technically, a finite small parameter e^ is introduced and the function is evaluated at a 
sequence of points tj = (j — ^)eh, then the coefficients are approximated to order 0(e^) 
by optimized finite differences. Obviously, parallel computation can be used to speedup 
the numerical e-expansion of hypergeometric functions. 

In the analytic e-expansion method, the coefficients of different orders are calculated 
separately. Usually, the expansions of A(e) and B(e) are needed to obtain the expansion 
of A{e)B(e). Somewhat differently, in the numerical e-cxpansion method, all coefficients 
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of the expansion arc extracted simultaneously. In addition, A(e)B(e) is expanded directly 
and there is no need to calculate the expansions of A(e) and B(e) beforehand. In other 
words, the numerical e-expansion method makes the evaluation of complex expression 
fast and convenient. 

The implementations of the algorithm arc presented in the form of a Mathcmatica 
package NumExp.m and a Python package numexp.py. Note that the mpmatri [37] package 
is required by the Python interface. The algorithm may also be implemented in other 
computation systems. Recently the NumExp package has been applied successfully to 
the doubly heavy hadron spectral density calculation [38], where q +iF q and/or Appcll 
function F4 are involved. 

The paper is organized as follows. In the next section the theoretical background and 
the design of the program are described. Then the usages of the packages are shown in 
Sec. 3. Whereafter, some examples are presented in Sec. 4. Finally, a summary is given 
in Sec. 5. 

2. Theoretical background 

2.1. Laurent expansion 

Formally, the Laurent expansion of the e-regularized hypergeometric function can be 
expressed as 

00 

f(z,e) = Y^fn(z)e n > keZ, (1) 

n—k 

where f(z, e) is a short form of the hypergeometric function, z the argument vector, and 
e the regularization parameter, respectively, e can be positive, negative or even complex. 
Practically, the new function F(z,e) = e~ k f(z,e) instead of f(z,e) is used to perform 
the Laurent expansion. If k < 0, this procedure makes the new function finite for e — > 0. 
On the other hand, if k > 0, this procedure helps to increase the precision of the results. 
So, the following Laurent expansion 

— 1 00 

F(z,e)= Y, e n -0 + Y,Fn(z)e n (2) 

n— — m n— 

will be used to derive the formulas of numerical epsilon expansion of hypergeometric 
functions. Note that the spuriously divergent terms y~] n —_ m e™ ■ vanish for analytic 
calculation, but multi-precision computation is needed to get rid of round-off errors from 
such terms for numerical calculation. 

The coefficients F n (z) are related to the partial differentials of F(z, e) with respect to 
e. In this work, the finite differences of F(z, e) will be used for the numerical calculation. 

2.2. Finite difference method 

When the regularization parameter e is set to zero, function F{z, e) is usually ill- 
defined. With non-positive integer parameter, Gamma function and hypergeometric 
function contain singularities and may not give desired results. Moreover, not only the 
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coefficient at the leading order in e is needed to be extracted. Therefore, some difference 
formulas like 



F (z) = F(z,e) + O(e), 

Fi{z) = F( Z ,2e)-F( Z ,e) +0(e)j 

Mz) = F( Z ,3e)-2F^2e) + F( Z ,e) + ^ 

and so on can be used to extract the coefficients of the Laurent expansion. That is, 
F(z, 0) will not be used in the calculation. 

The precision of the formulas above is very low. Actually, there are many ways to 
improve the precision and keep the computational complexity unchanged [36]. In this 
work, the algorithm 

n Qn 

w = E ~i F (*> I) + <^r +1 ) (3) 

3=0 

is used for the numerical epsilon expansion. C™- is the nth degree weight array of the 
finite difference method, and < i,j < n. It is trivial to work out these constant arrays 
by solving some linear equations. For example, 

3 3 
-1 1 

1 -2 

The computation of F(z, (j — h)^h) is t ne most time-consuming part of the whole 
numerical calculation. It is worth noting that only n + 1 times of function evaluation is 
needed to obtain all Fi( Z ) (i = 0, . . . , n) to order 0(e^~ l ). In particular, with 3 times of 
function evaluation, one obtains F (z) to order 0(e 2 h ), Fi(z) to order 0(e{) and F 2 {z) 
to order C(e°), respectively. If eh is numerically small enough and/or n is large enough, 
one could get the numerical e-expansion of F(z, e) to the desired precision. 



2. 3. Precision 

However, eh cannot be too small. From Eq. (2) and Eq. (3), it is easy to see that 
— (m+n+1) lg(\eh\) digits of working precision is needed to obtain Fi(z) to order 0(e h l ~ 1 ). 
If e/i is too small and the working precision is insufficient, the precision of the low order 
coefficients is restricted by the working precision and the high order coefficients will be 
inaccurate or meaningless. On the other hand, higher precision means lower speed of 
computation. Therefore, one has to find a balance between precision and efficiency. 

With a given working precision, one can use either smaller (and less function 
evaluations) or larger e% (and more function evaluations) to obtain the coefficients of the 
numerical e-expansion. For example, with 15 digits of working precision and Ch = 10~ 5 , 
Fq(z), Fi( Z ) and F 2 (z) have roughly 15, 10 and 5 digits of precision, respectively. If 
e h = lfT 3 , F (z), Fi(z), F 2 {z), F 3 (z) and F 4 (z) have roughly 15, 12, 9, 6 and 3 digits 
of precision, respectively. It is quite obvious that lower coefficients can be evaluated to 
high precision even with a large at the cost of more function evaluations. It also shows 
that the precision of the coefficients can be increased in a trivial way. 
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Multi-precision computation is the key to the success of the numerical e-expansion 
method. Mathematica and the Python library mpmath [37] provide algorithms to eval- 
uate hypergeometric functions p F q and lots of special functions to arbitrary numerical 
precision. In particular mpmath provides a function hyper2d ( ) to evaluate a wide variety 
of hypergeometric functions of two variables, such as Appcll functions Fx, F2, F3 and F4, 
which occur in the calculations involving massive particles [2, 39, 40, 38]. With the help 
of these computation systems, the package NumExp can be used to perform e-expansion 
of hypergeometric functions in a numerical and efficient way. 

3. Usage 

The algorithm in the preceding section is implemented in the package NumExp, and 
the package can be obtained from the ancillary files of arXiv: 1209.3971. The package 
provides both Mathematica and Python interfaces for users. 

3.1. Mathematica interface 

After a successful installation, the Mathematica package NumExp. m may be loaded 
with the command 

In[l] := « NumExp 1 

The package provides three public functions: 

• NumExpFDC[n] generates the nth degree weight array C™-. If n — 2, the array in 
Eq. (4) will be generated. In most cases, this function will not be used directly. 

• NumExp [expr , ep , eh, n] is the main function of this package. It expands expr to 
order C(ep n ). eh is the numerical step parameter in Eq. (3). The usage of this 
function is illustrated by the following example: 

In [2] : = f [e_ ,z_] : =Hypergeometric2Fl [1 ,-e , 1-e ,z] ; 

In [3] := NumExp[f [e,3] ,e,10~(-4) ,4] 

Out [3] = {{1 . 0000000000000000001 ' 20 , 1 20*1} , 

■CO . 693147180559944 ' 16 , 3 . 141592653589793 1 16*1} , 

{-2 . 32018042335 1 12 , 3 . 45139229523 ' 12*1} , 

-C-3.7421220'8, 1 . 8958710 ' 8*1} , 

{-3.7511'4, 0.6944'4*I}} 

The precision of the output is estimated by Mathematica automatically. Note that 
the precision of the high order coefficients is slightly overestimated, which is the 
effect of the 0(e^ +1 ) remainder terms. One can use smaller e/j or take higher order 
expansion to increase the precision of the coefficients. 

Besides, the precision of the input parameters must not be lower than the working 
precision of NumExp, or the precision of the result will be limited. 

• ListNumExp [1st , eh,n] extracts the coefficients of the e-expansion from a list of 
values /(eo), /(ei), . . ., f(e m ) which are evaluated beforehand, eh is the numerical 
step parameter and n < m. The code 
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In [4] := lst=Table[f [(j-l/2)*lCT(-4) ,3] ,{j ,0,6}] ; 
In [5]:= ListNumExp[lst,10~(-4) ,4] 



will generate the same results. If n < m, the first n + 1 (from 1 to n + 1) values 
of 1st will be used by ListNumExp. Note that the step parameter eh must be 
consistent with the parameter in the f(ej) evaluations, or the expansion will be 
nonsense. 

This function is especially useful because the time-consuming /(e 3 ) evaluations can 
be performed independently, which makes the parallel computation possible. 

3.2. Python interface 

The Python package numexp . py is specially designed for numerical e-expansion of 
hypcrgcomctric functions. Since e is used to regularize the expression, these functions 
may contain spurious divergences at e — and should not be evaluated at e = directly. 
In the package, the function is evaluated at a list of points (j — \)^h and the coefficients 
are evaluated by finite differences of these values. 

It is worth noting that mpmath [37] provides a function taylor () to produce a degree- 
n Taylor polynomial around the point x of the given function f(x). The differentials are 
approximated by finite differences and this function can handle singularity of /(0) by 
shifting the points half of a step length h/2. However, taylor () is pretty slow when the 
precision and the order of expansion are high because the default step parameter h is 
related to the working precision and is too small. If a fairly large h is specified, the speed 
increases while the precision decreases. Moreover, if the singular option is set to True, 
which is needed by the expansion of regularized hypergeometric functions, order 0(n 2 ) 
instead of order O(n) algorithm is used by taylorO and the calculation is roughly n 
times slower. 

numexp. py always uses 0(n) algorithm to calculate the coefficients of the expansion. 
With the given number of function evaluations, every coefficient is calculated to the best 
precision. The error is controlled in a systematic way and the precision is guaranteed as 
mentioned in Sec. 2.3. Therefore, numexp. py is recommended for numerical e-expansion 
of the regularized hypergeometric functions. 

numexp. py provides two public functions: 

• numexp (f , eh, n, args= [] ) expands f (ep,*args) to order C(ep n ). args is the 
argument list of function / and eh is the numerical step parameter. The usage of 
this function is illustrated by the following example: 

In[l]: from numexp import * 
In [2] : mp . dps = 20 

In [3]: f2 = lambda ep: hyp2f 1 (1 , -ep, l-ep,3) 
In [4]: numexp(f 2,mpf ('le-4') ,4) 
Out [4] : 

[ (1 . 0000000000000000001-6 . 2579670696103073837e-21 j ) , 
(0 . 69314718055994379081+3 . 1415926535897933231 j ) , 
(-2.320180423351603576 +3 . 451392295225348322 j ) , 
(-3 . 7421219941535409023+1 . 8958709608982570633j ) , 
(-3. 7510818776117221586+0. 69441875410183733333j )] 
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These five coefficients have roughly 20, 15, 12, 7 and 3 digits of precisions. One 
can take the first n — 1 coefficients and discard the last one or two. This may be 
an appropriate strategy for the e-cxpansion. 

If the function /() contains more than one argument, all other arguments except e 
should be put in a list. 

In[5]: f3 = lambda ep,a,z: hyp2f l(a,-ep, l-ep,z) 
In [6] : numexp(f3,mpf ('le-4') ,4, [1,3]) 

• listmimexp(lst,eh,n) extracts the coefficients of the e-expansion from a list of 
values /(eo), /(ei), . . ., /(e m ) which are evaluated beforehand, eh is the numerical 
step parameter and n < m. If n < m, the first n + 1 (from to n) values of 1st 
will be used by default. 

In the following of this work, only accurate digits of the result will be kept in the 
results for the sake of simplicity. 

4. Examples 

In this section some examples for the numerical e-expansion of hypergeometric func- 
tions are given. So as to keep the expressions short, fairly large e/j is used in the follow 
calculation. If the coefficients of the e-expansion are needed to be evaluated to a higher 
precision, smaller e/j should be used. 

4-1- Hypergeometric functions and analytic continuation 
The simple Gauss hypergeometric function 



can be analytically expanded by several packages [26, 28]. Alternatively, it is easy to 
expand this function by NumExp in a numerical way. For z = 0.3, one gets 



A(0.3, e) =0.999999999999999999 + 0.326129510 e 2 + 0.737565 e 3 + C(e 4 ), (6) 



where the input parameters e/j = 10~ 4 and n = 4 are used and only accurate digits are 
kept in the result. 

If the hypergeometric function is divergent for e — > 0, negative power terms of e occur 
in the Laurent expansion. For the function 



one can use eB(z 7 e) instead of B(z,e) to extract the coefficients, as stated in Sec. 2.1. 




(5) 




(7) 



Explicitly, 



5(0.5, e) =4.00000000000000000000 e" 1 4- 0.9999999999999999 
+ 6.408403647539 e - 10.5952025 e 2 + C(e 3 ). 
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(8) 



If \z\ > 1, analytic continuation is used to evaluate the hypergeometric function. 
Numerically, one gets 



,4(3, e) =0.999999999999999999999+ (2.320180423313- 3.451392295223 i)e 2 
- (5.64797470 + 12.27713795 i)e 3 - (30.5240 + 8.050 i)e 4 + C(e 5 ), 



(9) 



where analytic continuation is performed automatically in the numerical calculation. 

It is worth noting that analytic continuation is a non-trivial issue for the analytic 
expressions generated by HypExp and others [41, 29]. If one directly set z = 3 in the 
analytic expansion of Eq. (5), the coefficient of e 4 term will be —30.5241 — 7.35203 i, 
where the imaginary part is incorrect. Technically, a small imaginary part should be 
introduced into the argument to keep it on the correct side of the branch cut. That is, 
one can use z = 3 — 10~ 20 i to reproduce the correct coefficient. 

There is another issue on the analytic continuation of the hypergeometric functions. 
For the q +\F q functions, the analytic continuation formula for \z\ > 1 is 



q+lF q 



ax, . 



7 a q+l 

■ ■ , b„ 



r(ai)---r(o 9+ i) 



E 



X (-*)" 



q+ 



n?=irfe-0 

a>i, {1 + a-i - bk}k=i,...,q 

{1 + a, — ak}k=l,...,q+l;kjti 



, (10) 



where 



i z. if 



(ii G Z, one can introduce another auxiliary regularization 
parameters e' into the parameters a.^s to regularize the analytically continued expression. 
With a,i — > a.i + ait' and o.j — on ^ 0, the analytically continued expression is well 
regularized. Technically, it is advisable to use irrational ctj's to fulfill the requirement of 
aj — OLi 7^ 0. Note that the expression is finite in e' and e' can be set to a small enough 
value in the numerical calculation. 

In Rcf. [38] the doubly heavy hadron spectral density has been expressed in the form 
of e-regularizcd hypergeometric functions. A part of the correlation function is 



nf rt (z) 



3 f r(-e-4)r(-e- 2) 2 r(-e) 
647T 6 \ r(2-e)r(-2e-4) 

'-6-4,-6-2 



X 3^ 



2 2 

9 ) ^ 



r( - e -4)r(-6-2) 2 r(-6) 



x 3 F 2 



r(3-e)r(-2e-4) 
'-e-4,-e-2 



1,3 



V2e',-e 
r 



(11) 



The auxiliary e' is introduced to regularize the analytically continued expression. In 
Rcf. [10], a similar regulator 5 was used to cancel out the spurious divergences. By using 



Pi?) = — Im{n(^ + ie) 



II(z — is)} = ImII(z — ie) 

TT 



(12) 



and setting e' = 10~ 10 , it is easy to extract the spectral density from the e-expansion 
of Eq. (11). Although the Laurent expansion of the correlation function n(z) contains 
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e~ m terms, their coefficients are all real and the imaginary parts are zero. Actually, one 
can directly evaluate the imaginary part of Eq. (11) to order 0(e ) and no numerical 
e-expansion is needed at all. Even so, one can use NumExp to improve the precision of 
the result. 

Another unorthodox regularization instead of the dimensional regularization has also 
been used in Rcf. [38] to obtain well regularized hypergeometric representation of the 
spectral density p(z), where all coefficients of e~ m terms are exactly zero. Then one can 
use a small e to evaluate the special density directly. 

4-2. Generalized hypergeometric functions 

Integral representations of some hypergeometric functions may be used to perform the 
e-expansion. In Ref. [42], the program SecDec has been used to evaluate the e-expansion 
of two generalized hypergeometric functions 



A{e) and -B(e) can be analytically expanded by HypExp [28, 29]. SecDec [42] takes several 
minutes to obtain the coefficients. It is transparent that multi-dimensional numerical 
integration is time-consuming and its precision is low. 

The numerical epsilon expansion method might also be applied to the e-regularizcd 
parameter integrals, if the poles are subtracted properly. However, the precision of the 
multi-dimensional numerical integration is low, which makes it difficult to extract high 
order coefficients of the epsilon expansion by the finite difference method. 

The hypergeometric representation has the advantage in high precision calculation. 
NumExp takes less than one second to expand this type of generalized hypergeometric 
functions to order C(e 6 ). The numerical results are 



where e/j = 10~ 4 and n — 6 are used to perform the numerical e-expansion. Note that 
all inaccurate digits of the coefficients are discarded. The precision of the results can be 
improved by using smaller e^ and/or larger n. 

The analytic e-expansion of generalized hypergeometric functions about rational pa- 
rameters was studied in Rcf. [24]. So as to illustrate the capability of NumExp, the numer- 
ical e-expansion of a hypergeometric function with rational and/or irrational parameters 
is performed here. For example, the function 




(13) 



(14) 



A(e) =0.9999999999999999999+ 0.189532432184360 e 

- 2.2990427423 e 2 + 55.469019 e 3 - 1014.39 e 4 + C(e 5 ), 
B(e) =0.999999999999999999- 4.27968776167885 e 

- 26.6975474079 e 2 + 195.87119 e 3 - 7313.7 e 4 + C(e 5 ), 



(15) 



(16) 




(17) 
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cannot be expanded by HypExp or other analytic e-expansion algorithms, while SecDec 
may be able to expand this function by using its integral representation. It is not sur- 
prising that NumExp can handle this kind of function with no difficulty. Explicitly, 



C(0.5,e) =1.0000000000000000000- 1.44555526747927 e 

+ 3.9383879447 e 2 - 266.94735 e 3 + 298.66 e 4 + C(e 5 ) 



(18) 



where e% — 10 4 and n — 6 are used in the expansion. 

4-3. Hypergeometric functions of two variables 

Hypergeometric functions of two variables occur in the Feynman diagram related 
calculations of three kinematic variables and/or masses [2, 39, 43, 40, 38]. Recently, a 
few algorithms or packages [25, 26, 27, 31] have been developed, which can be used to 
perform the analytic e-expansion of some special types of hypergeometric functions of 
two variables. To the authors' knowledge, these implementations can only be applied to 
very limited cases. 

The Python library mpmath [37] can evaluate any of the 34 distinct convergent second- 
order Horn-type hypergeometric series. With the help of mpmath, NumExp can be used to 
perform the numerical e-expansion of a wide variety of hypergeometric functions of two 
variables. 

The analytic e-cxpansion of the Appcll function 



is presented in Ref. [25, 26]. Here NumExp is used to show the possibility of the numerical 
e-expansion of such kind of functions. Explicitly, 



A(0.3, 3.4, e) =1.428571428571428571 - (2.43800237269678 + 4.48798950512827z)e 

- (7.0450867025 - 6.6246225829 i)e 2 + (6.57967 + 6.095515 i)e 3 , (21) 

where e/j = 10~ 4 and n = 4 are used in the expansion. For |x| + \y\ > 1, analytic 
continuation is performed automatically in the numerical evaluation. The parameters of 
this Appcll function above are pretty simple and the expression of the analytic expansion 
is short. If the parameters arc quite complex, the expression will become lengthy and 
the numerical e-expansion may be a better choice for the practical calculation. 

Presently, the analytic e-expansion of Appell function F4 is valid only for some specific 
cases. In Ref. [10], F4 was converted to 2-P1 or Fi, and then XSummer [27] was used to 
perform the e-expansion. Practically, Appell F4 and other hypergeometric functions of 
two variables can also be expanded by NumExp. 

In the doubly heavy hadron spectral density calculation [38] , Appcll F4 occurs if two 
heavy quarks arc different. For the (Qq)o(Q'q)o molecular state in Ref. [44], a part of 




(19) 



A(0.3, 0.4, e) =1.428571428571428571 + 0.700889880640673 e 
+ 1.6060080586 e 2 + 1.418379 e 3 + C(e 4 ), 



(20) 
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the correlation function can be expressed as 

np "<" r < 2 « - 4 > r <- e - 2 > r < 2 - e » f > (* ; - rr 2 



9m?- £ 7?4 +e , N , , /e-2;-2e 
2 -r(e - 2) 2 r(-2e)F 4 ' 



256tt 6 v ; v ; V 3 - e, 2 



s 

ml 



(22) 



where the Appell F4 functions are regularized in an unorthodox way. 

For m\ = 1, m 2 = 2 and s = 19, for example, the prefactors of the first and the 
second terms in Eq. (22) can be expanded as 

Boi(e) = - 1(T 5 x + + 81.7264 + 76.8361 e) + 0(e 2 ), (23) 

r, / s ^5 /7.31364 27.0103 95.1100 „ , r \ m . , 

B 02 (e) = - 10- 5 x f — — + —- 2 — + + 251.522 + 0(e), (24) 

respectively. In the same way, the Appell F4 functions of the first and the second terms 
in Eq. (22) can be expanded as 

Bii(e) = - 0.75 <T X + 47.3086 + (2.27125 - 6.55331 i)e 

- (93.1352 - 8.50940 i)e 2 + 0(e 3 ), (25) 

B l2 {e) =1.0 + 1.51042 e + (17.8362 + 8.73770 i)e 2 

+ (0.87460 - 34.6270 i)e 3 + 0(e 4 ), (26) 

respectively. -Boi(e) need to be expanded to e 1 because -Bn(e) contains e" 1 contribu- 
tion. Note that analytic continuation of F4 is performed automatically in the numerical 
calculation. Then, 



lP crt (19) = £ i(e)Bii(e) + B m {e)B l2 {e) 
0.00492267 0.00676321 



(0.0380542 + 0.0000352703i)- (27) 



It is transparent that the imaginary parts of e~ m terms are all zero. In fact, one can 
directly evaluate the imaginary part of Eq. (22) to order 0(e°) and no numerical e- 
expansion is needed at all. Even so, NumExp can be used to improve the precision of the 
result. Moreover, the spectral density p(z) can also be expressed by Appell F 4 functions 
and these e-regularizcd functions can be evaluated numerically. 

5. Summary 

The Feynman integrals can be calculated in a less tricky and systematic way by 
using the hypergeometric representation. Presently, the analytic e-expansion methods 
of hypergeometric functions have been used in practice. However, the application of 
such expansion has lots of limitations. For example, only some classes of hypergeometric 
functions with the specific form of parameters can be expanded, and some newly defined 
special functions are used as primary elements and the expanded expressions are quite 
lengthy which make the physical analysis cumbersome. 
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In this work it is demonstrated that the well-regularized hypergeometric functions 
can be evaluated directly and numerically. An algorithm is developed and the package 
NumExp is presented for expanding hypergeometric functions and/or other transcendental 
functions in a small rcgularization parameter. Hypergeometric function is expressed 
as a Laurent series in the rcgularization parameter and the coefficients arc evaluated 
numerically by using multi-precision finite difference method. This elaborate expansion 
method works for a wide variety of hypergeometric functions, such as the generalized 
hypergeometric functions p F q , the Appell hypergeometric functions and the Horn- type 
hypergeometric functions of two variables, which are needed in the context of dimensional 
regularization for loop integrals. The divergent and finite parts can be extracted from 
the final result easily and simultaneously. In addition, there is almost no restriction on 
the parameters of hypergeometric functions. 

The numerical e-cxpansion method may not be suitable for the intermediate expres- 
sions, but it is good enough for the final results. Practically, the divergent and finite 
parts can be extracted from the final result easily and simultaneously. Moreover, par- 
allel computation can be used to speedup the numerical e-expansion of hypergeometric 
functions. In other words, the numerical e-expansion method makes the evaluation of 
Feynman integrals fast and convenient. 

It is worth noting that this method can serve to be an important cross-check of the 
analytic e-expansion methods and it is possible to extend the method to expand functions 
with two or more regularization parameters. 
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